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Alwtrset. A one«parsineter family of explicit and implicit second-order* 
accurate, entropy satisfying, total variation diminishing (TVD) schemes 
has been developed by Harten. These TVD schemes have the property of 
not generating spurious oscillations for one-dimensional nonlinear scalar 
hyperbolic conservation laws and constant coefficient hyperbolic systems. 
Application of these methods to one- and two-dimensional fluid flows con- 
taining shocks (in Cartesian coordinates) yields highly accurate nonoscil- 
latory numerical solutions. The goal of this work is to extend these methods 
to the multidimensional Euler equations in generalised coordinate systems. 
Some numerical results of shock waves impinging on cylindrical bodies are 
compared with MacCormack’s method. 


§1. Motivation and Objeetivo 

Several techniques for the construction of nonlinear, second-order- 
accurate, high-resolution, entropy satisfying schemes for hyperbolic conser- 
vation laws have been developed in recent years. See, for example, van Leer 
[1], Colella and Woodward [2], Harten [3,4], Roe [5] and Osher [6]. We can 
also view these schemes as shock-capturing algorithms based on either an 
exact or approximate Riemann solver. From the standpoint of numerical 
analysis, these schemes are TVD for one-dimensional nonlinear scalar hyper- 
bolic conservation laws and for one-dimensional constant coefficient b^er- 
boiic systems. In [4], Harten introduced the notion of TVD schemes. Entropy 
satisfying TVD schemes have the property that they do not generate spurious 
oscillations and that the weak solutions are physical ones. The goal of con- 
structing these highly nonlinear schemes is to simulate complex flow flelds 
more accurately (i.e., to construct schemes that are stable in a strong non- 
linear sense). TVD schemes are usually rather complicated to use compared 
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with the conventional shock-capturing methods such as variants of the Lax- 
Wendroflf scheme. The complexity of these schemes has inhibited their ap- 
plication to complicated flow geometries in the past. 

Application of Harten’s explicit and implicit methods to standard one- 
and two-dimensional transient and steady-state gas-dynamic test problems 
in Cartesian coordinates was examined by Harten (3] and Yee et al. [7-9J. 

In both one and two dimensions, accurate solutions containing shocks and 
contact discontinuities were obtained. 

The objective of this report is to extend Harten’s TVD method to general- 
ized coordinate systems, and to test the method on a two-dimensional 
problem of a moving shock wave impinging on a cylinder. The numerical 
results are compared with MacCormack’s explicit method. From here on, we 
refer to this method as the TVD scheme. 

A description of the TVD algorithms in Cartesian coordinates can be found 
in reference [9]. A description of this method for two-dimensional Euler 
equations in generalized coordinate systems will be discussed in the next 

section. Some results on the shock wave-cylinder interaction are presented in 
section 3. 


§2. Extension of an ExpUcit TVD Scheme for the Euler Equations in 
Genefaliied Coordinate Systems 

Here we assume the reader is either familar with the development and 
propeHies of explicit and implicit TVD schemes, or will consult the ap- 
propriate references |3,4,7-9]. A brief description of these methods and a 
detailed implementation of these methods for one- and two-dimensional Euler 
equations of gas dynamics can be found in reference [9]. To avoid extra 
notation, a particular form of the explicit TVD scheme in [3] is extended 
to generalized coordinates. Generalization of the implicit TVD scheme to 
arbitrary geometries follows the same procedure. 


S2.1 The Euler Equations 

% 

In two spatial dimensions, the Euler equations of gas dynamics can be 
Written in the conservative form as 
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with m — pu and n = pv. The primitive variables are the density p, the 
velocity components u and v, and the pressure p. The total energy per unit 
volume e, is related to p by the equation of state for a perfect gas 


P = (7“l) 
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(m^ w^) 

2p 
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where 7 is the ratio of specific heats. 

A ^neralized coor^nate transformation of the form ^ = ^(ar,y) and rj = 

^hich maintdns the strong conservation law form of equation ( 1 ) is 
given by 


dQ 9f($) ag($) 

d( dti 


( 2 ) 


where 0 = Q/J. f" = {(.F + 6 = + „,G)/J, and / = 

— ^yVx, the Jacobian of transformation. Let A = dF/dQ and B = 
BG/dQ; then the Jacobian A and ^ of and d can be written as 




Let c be the iocal speed of sound; the eigenvaiues of A are 
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where tj = + ^yv) and + e®. The eigenvalues of B are 

(<»44<) = (^“ y, V-^k^c, V) (4b) 

where V — -f rtyv) and ^ = ^2 _|_ ^2 

Furthermore, let = {R\,R%R%R\) and R^ = {R\,R%,R^^,R*) be 
the matrices whose columns are eigenvectors of A and B. Let R^^ and R ~^ 
be the inverses of R^ and R^f. A form of R^ and can be written as ’ 
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Let the grid spacii^ be denoted by and A;? such that ( = and 
*7 = kAr}. Denote Oj+1/2,* as some symmetric average of < 5 ^,* and Qj+i,k 
(for example, the Roe’s average [10]). Let aJ_|.i/2, ^j+1/2, ^J+1/2 <Jenote 
the quantitic'^ of a^, R^, related to A evaluated at <^>+1/2,*- Similarly, 
let ^*+1/2, Rr+i/2 denote the quantities of aj,, R^, R~^ related to 

B evaluated at ^j,k+i/2^ 

We define 


ay+i/2 = ^j+i/2(Qi+ i.fc - Q},k) (7 a) 


A A 

as the component of (Qy+i,* — Qj,k) (omitting the k index) in the locally 
/•th characteristic f-direction [ 9 ]. Denote 


«*+l/2 = Rk-ll/2iQhk+l - Qj\k) 


( 7 b) 


as the component of {Qj,k+i — Qj,k) (omitting the j index) in the locally 
/•th characteristic 17-direction. The vector a of equation ( 7 a) can be written 

as 
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where 
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cc = — AriA,+i/2n — (A;2W/+i/2 — feit^i+i/2)Aj+i/2/> + k2Aj+i/2m (8d) 
irith 

^j+i/2p — Pj+i,k — Pj,k, Aj+i/2m =s mj4-i,* ~ mj,* (8e) 

and 

Aj+i/2n = nj+i,* — n>, *, Aj+i/2C = c^+i,fc — gj,* (8f) 

The simplest form for Qj+i/2,k is 

$i+i/2,fc = (<5/+i,fc + 0j,*)/2 (9) 

Roe’s form of the averaging in the (-direction is: 
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Therefore to use Hoe*s averaging for the ^-differencing, ali one has to do is 
compute Uj+i/ 2 ,k, Vj+i/a,*, and cj+i/ 2 ,k in equations (4a), (5)-(6), (7a) and 
(8) by equation (10). Similarly, Roe’s averaging can be obtained for Uj^k -^i/ 2 > 
* 4 - 1 / 2 ) nnd Cj,k+i/2> in the numerical experiments for the two-dimensional 
test problem. Roe’s averaging is used. 
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§2.2 Algorithm of a TVD Method In Ueneraliaed Coordinates 

i ® particular form of the explicit TVD 

method of (3J in generalized coordinates, when implemented by the method 
of fractional steps, can be written as follows: 
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We can define the variables of equation (lie) by simply replacing the 
subscripts J and / + 1/2 in equations (llf)-(l ll) by A: and A; + 1/2. Extension 
of the implicit TVD scheme [4,9] to a generalized coordinate system follows 
the same procedure as described above. 

In general, if one handles the intermediate boundary conditions correctly, 
one only needs to do the half steps in (11a) at the beginning and immediately 
before printout; i.e., 







where 



is the initial condition. 
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$3. Numerical Result for the Shock Wave-Cylinder Interaction 

A good test problem for assessing the capabilities of any shock-capturing 
scheme Is the shock-diffi'action problem; i.e., the computation of the unsteady 
flow field resulting from a planar moving shock wave striking an obstacle. In 
the present numerical experiment the difilraction process is determined over 


a cylinder. The shock patterns at two Instances in time and after initial 
ifiipingement are sketched in Fig. 1. 

When the Incident shock first collides with the cylinder, regular reflection 
occurs at the shock impingement point. As the impingement point of the 
Incident shock propagates around the body, the reflection process makes a 
transition from regular to Mach reflection. It rhould be pointed out that 
during the transition process, complex and double Mach reflection shock 
structures are possible. Their occurrence is dependent on the initial strength 
of the incident shock wave. For single Mach reflection, a triple point forms 
and the incident shock no longer touches the body. Emanating from the triple 
point are three waves: 1) a Mach stem which strikes the body perpendicularly, 
2) a slip surface or shear layer Which strikes the body and results in a 
vortical singularity (nodal point of streamlines), and 3) the reflected shock 
which propagates away from the body. In addition to the above flov field 
characteristics, a stagnation point (saddle point of streamlines) exist.'. *.he 
plane of symmetry, both forward and aft on the body. 

The shock-diffraction problem contains most of the fl; v i.scoiiitl*:. , ilies pos- 
sible with the Euler equations and is thus a test or a numerical shock- 
capturing procedure. Both MacCormack’s e>piicit method and the explicit 
TVD scheme were applied to the shock wave-cylinder interaction problem. 
For afair comparision, the TVD scheme was implemented in an existing com- 
puter code [11] which also contained MacCormack’s method, so that same 
initial conditions, boundary conditions and coordinate transformation were 

tt^. A cylindrical grid consisting of 50 points around the half-cylindrical 
(^-direction) boojr and 51 points between the body and outer boundary (in- 
direction) was 'ised. The body radius is one and the distance from the body 
to the outer ^undary is between 2 to 4 (depending on the incident shock 
Mach number). Rays from the coordinate system origin are spaced at equal 
angles with points uniformly placed in the radial direction between the body 
and the outer boundary (see Fig. 2). 


§3.1 Initial and Boundary' Conditions 

Fig. 2 shows a schematic representation of the grid with its boundaries 
and initie! conditions. The nodal points to the right of the planar moving 
shock are Initialized to free stream values while those to the left are set equal 
to the post moving shock conditions. In the outer boundary, it is necessary 
to track the moving planar shock as a function of time along this boundary 


surface. 


At the planes of symmetry, the reflection principle is used; i.e,, the pressure, 
dciiLr’.'^ niid u-veiocity component are treated as even functions across the 
plane of symmetry while the u-velocity component is treated as an odd 
function. The boundary condition at the surface of the cylinder must satisfy 
the tangency condition which requires that the velocity in the radial-direction 
be equal to zero at the body. Furthermore, for convenience, an image line of 
nodal points is considered which falls one mesh interval inside the body, so 
that-the reflection principle can be applied. 


§3.2 Numerical Result 

MacCormack*s method with a fourth-order dissipation term was run at 
a Courant (CFL) number of 0.6 for stability while the TVD method was 
operated at a Courant number of 0.9 for efficiency. The Courant number 
is a measure of the maximum permissible time step for a stable solution. 
The TVD method is insensitive to Courant number between 0.5 and 1. Both 
methods were run to approximately the same total time (100 steps for the 
MacCormack’s method, 70 steps for the TVD method). The results in the 
form of pressure and density contour plots are shown in Fig. 3 at a time for 
which Mach reflection of the incident shock exists. The incident shock Mach 
number was 2. The results from MacCormack’s method are shown in Figs. 
3a and 3b. Those for the TVD method are shown in Figs. 3c and 3d. It 
can be seen that the IVD scheme results in a better defined flow field; i.e., 
“crisper” shocks and hardly any associated spurious oscillations. The slip 
surface which emanates from the triple point is smeared beyond recognition 
by both methods. It is, however, possible to observe the location on the body 
where the slip surface impinges (i.e., the vortical singularity). At this point, a 
local pressure minimum occurs, and the pressure contours, as a result, encircle 
it. This behavior can be observed in the pressure contour plot in the region 
to the left of the Mach stem at the body for the TVD scheme. 

To test the shock-capturing capability of the TVD method at higher inci- 
dent shock Mach numbers using the same grid size as before, cases were run 
at Mach numbers between 3 and 10. Figure 4 shows the pressure and density 
contour plots for the TVD scheme with an incident shock Mach number of 
10 at a time when the incident shock was already passed the cylinder. For 
this case the Mach stem that extended from the triple point to the body 


hM .truck the plune of symmetry and reUected from it. This will eveatually 
™ M It did at the body. The result 
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§4 Concluding Remarks 

“woo^order accurate explicit TVD scheme in generalieed 
coordinate systems has been applied to obtain transient soiutions on the two- 

“ moving shock wave impinging on a cyiinder. Fairiy 
*®" obtained. Moreover, from numericai experiments 
k 's stabie in a strong nonlinear sense (e.g., the calcuiation with an 

the TTO^brw^^ “'"“bet of 10). The report is the first attempt to apply 
K* ““-'“ttesian coordinates. It is preliminary in natS^ 

** underway on improving the resolution of discontinuities bv 

5m^U?te„r‘oTM" ■“!'«*'>'> t‘2-131, bnd on improving the eSc“ o^f 
computation by possibly using a large-time step (explicit method) approach 
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Fig. 1. Shock structure for shock diffraction 
over cylinder at tiivo different times. 
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Fig. 2. Schematic of the computational grid. 






cylinder interaction. {Ms = 10). 



